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We study the quantum phase transition in the two-dimensional random Ising model in a transverse 
field by Monte Carlo simulations. We find results similar to those known analytically in one- 
dimension. At the critical point, the dynamical exponent is infinite and the typical correlation 
function decays with a stretched exponential dependence on distance. Away from the critical point 
there are Griffiths-McCoy singularities, characterized by a single, continuously varying exponent, 
z' , which diverges at the critical point, as in one-dimension. Consequently, the zero temperature 
susceptibility diverges for a range of parameters about the transition. 
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Though classical phase transitions occurring at finite 
temperature are very well understood, our knowledge of 
quantum transitions at T — is relatively poor, at least 
for systems with quenched disorder. There is, however, 
considerable interest in these systems since they (i) ex- 
hibit newptmiversality classes, and (ii) display "Griffiths- 
McCoy"™^! singularities even away from the critical point, 
due to rare regions with stronger than average interac- 
tions. 

Just as the simplest model with a classical phase transi- 
tion is the Ising model, the simplest random model with 
a quantum transition is arguably the Ising model in a 
transverse field whose Hamiltonian is given by 
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Here the {erf} are Pauli spin matrices, and the nearest 
neighbor interactions and transverse fields hi are both 
independent random variables. This model should pr 
vide a reasonable description of the experimental systeml 
LiHo 3; Yi_ 2 ;F4 and may also be an appropriate modelQ to 
describe non-fermi liquid behavior in certain /-electron 
systems. 

Naturally the random transverse field Ising model has 
been quite extensively studied.-and many surprising an- 
alytical results are availabletTLl for the case of dimen- 
sion d = 1. For example, the dynamic critical expo- 
nent, z, is infinite. Instead of a characteristic time scale 
£ T varying as a power of a characteristic length scale £ 
according to £ T ~ £ z , one has instead an exponential 
relatione! £ T ~ exp(const. where ip — 1/2. This is 
called activated dynamical scaling. In addition, distri- 
butions of the equal-time o~f—o~f, r correlations, are very 
broad. As a result average and typicala correlations be- 
have rather differently, since the average is dominated 
by a few rare (and hence atypical) points. At the crit- 
ical point, for example, the average correlation function 



falls off with a power of the distance r as C av (r) ~ r , 
whereB fj = (3 — y5)/2 ~ 0.38 whereas the typical 
value falls off much faster, as a stretched exponential 
Ctyp(r) ~ exp(— const. r a ), with a = 1/2. As the criti- 
cal point is approached, the average and typical correla- 
tion lengths both diverge but with different exponentsB, 

1- e. £av ~ <5 _1/av ; £typ ~ <5~ i/typ , where S is the deviation 
from criticality, and v av = 2, vtyp = 1. Finally, there are 
strong Griffiths-McCoy singularities at low temperature 
even away from the critical point, coming from rare re- 
gions which are "locally in the wrong phase" . These are 
characterized by a single continuously varying exponentfl, 
z'(S), which diverges as S — > 0. 

An important question is whether these striking an- 
alytical results are a special feature of 1-d, or whether 
they are valid more generally. Unfortunately, the ana- 
lytical approach is only valid in 1-d, and very little is 
known in higher dimensions. Senthil and SachdevH9 have 
studied the model in Eq. (0) with site dilution and shown 
that activated dynamical scaling occurs along that part 
of the zero temperature phase boundary which is pre- 
cisely at the percolation concentration. However, it is 
not clear if this result also holds for the rest of the phase 
boundary, and, to our knowledge, there are no results 
at all for other, more general, models. Here, we inves- 
tigate the behavior of the random transverse field Ising 
ferromagnet in two dimensions by performing large-scale 
Monte Carlo simulations. Because the ferromagnet has 
no frustration we are able to use highly efficient clus- 
ter algorithms which considerably reduce critical slowing 
down. Our main conclusion is that the behavior of the 

2- d system is very similar to that of 1-d. 

In order to capture the random quantum critical be- 
havior in the intermediate size systems that we can sim- 
ulate, we wish the disorder to be effectively quite strong. 
In particular, we would like some of the fields to be much 
stronger than the bonds in their vicinity and vice-versa, 
which is captured by having distributions for both the 
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fields and interactions with a finite weight at the origin. 
We therefore take the following "box" distribution 



1, for < Jij < 1 
0, otherwise 

h^ 1 , for < hi < h 
0, otherwise 



(2) 



As is standardQ, we represent the d-dimensional quan- 
tum Hamiltonian in Eq. ([!]) by an effective classical ac- 
tion in (<i+l)-dirncnsions, where the extra dimension, 
imaginary time, is of size (3 = 1/T, and is divided up 
into L T = (3/ At intervals each of width At in the limit 
At — ► 0. The partition function can then be written 
as Z = liniAr^o Trexp(— S), where the effective classical 
action is given by 



£ K iS Si (t)S 3 (t) - ]T KiSi {r)St (t') 

(«,i),T i,r 



(3) 



where, r' = 



At, Si{r) = ±1, 
and exp(— IK/) 



tanh(AT^). (4) 



To study large systems sizes with small statistical er- 
rors, we use cluster algorithms which simultaneously flip 
many spins. For our results on Griffiths-McCoy singulafk 
ties we have developed^ a variant of the loop algorithir£3 
in which the required limit At — > is explicitly taken. 
We shall call this the continuous imaginary time algo- 
rithm. It represents the original quantum Hamiltonian 
exactly (apart from statistical errors). We tune through 
the transition by varying h. 

In our simulations which determine for critical expo- 
nents, we use a different approach, and exploit univer- 
sality according to which the universal quantities should 
be independent of At and so, for convenience, we set 
At = 1. We then have a three-dimensional Ising model, 
with disorder perfectly correlated in one direction, which 
we simulate using the Wolflo cluster algorithm. We shall 
call this the discrete imaginary time algorithm. It does 
not represent the quantum Hamiltonian exactly but is ex- 
pected to be in the same universality class. Furthermore, 
for this algorithm we find it convenient to parameterize 
the strength of fluctuations by an effective classical tem- 
perature T c i = 1 /Pel, (not equal to the real temperature 
which is the inverse of the size in the time direction) and 
write Z = Tr exp(— P c (H c i) where 

H cl = - Ji 3 S 1 (t)S j (t) - ^ JiSii^Siir + 1). (5) 
(ij),-r 1 > T 

Here t runs over integer values, 1 < r < L T , and the 
distributions of the interactions are given by 



1 for < Jij < 1, 
otherwise, 

2exp(-2J i ) for Jj > 0, 
for Ji < 0, 



which, from Eq. (0), is similar to Eq. (g). We tune 
through the transition by varying the classical temper- 
ature, T c i. 

The lattice is of size L in the space directions, and we 
denote disorder averages by [■ • -] av and Monte Carlo av- 
erages by (•••). For both algorithms we employ periodic 
boundary conditions in all directions. 

First of all we discuss our results for the location of 
the critical line using the continuous imaginary time al- 
gorithm. We do this by computing the Binder ratio 



,9a 



(M 4 

W: 



where M = J Q Si(r) dr, At any finite temperature 
the system is expected to be in the universality class of 
the classical two-dimensional classical random bond Ising 
ferromagnet. For small temperatures the size of the clas- 
sical critical region shrinks and we need to study larger 
sizes (we went up to L = 32) to get a reliable estimate of 
h c (T). By extrapolating the latter to T = 0, see Fig. [I], 
we obtain for the location of the quantum critical point 
h c = h c (T = 0) = 4.2 ±0.2. 
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FIG. 1. Results obtained with the continuous imaginary 
time algorithm. The left hand axis indicates the phase di- 
agram of the d = 2 random transverse Ising model: PM 
means paramagnetic, FM means ferromagnetic, T c — 1.00(1) 
is the critical temperature of the classical random Ising fer- 
romagnet (hi = 0) with the box bond distribution in Eq. (Q), 
and h c — 4.2(2) the location of the quantum critical point 
we are interested in. The right hand scale indicates the val- 
ues of d/z'(h) obtained from analyzing the integrated prob- 
ability distribution of lnxiocai according to (^) in the Grif- 
fiths-McCoy region, h > h c . The open circle corresponds to 
z'(h c ) = oo and the horizontal line at d/z' — 1 indicates the 
expected limit lim^^oo z'(h) — d. The broken line is just a 
guide to the eye. 

Now we turn our attention to the Griffiths-McCoy re- 
gion in the disordered phase (h > h c ). Due to the pres- 
ence of strongly coupled regions in the system the prob- 
ability distribution of excitation energies (essentially in- 
verse tunneling times for these ferromagnetically ordered 
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clusters) becomes extremely broad. As a consequence we 
expect the probability distribution ofi-local susceptibili- 
ties to have an algebraic tail at T = OE^rta 

^(ln Xlocal) ~ 7777 In Xlocal, (8) 

z '(h) 

where f2(lnxiocai) is the probability for the logarithm 
of the local susceptibility \i a * site i to be larger than 
In Xiocal- The dynamical exponentu z'{h) varies continu- 
ously with the distance from the critical point and param- 
eterizes the strengths of the Griffiths-McCoy singularities 
also present in other observables. At finite temperatures 
the distribution of Xlocal is chopped off at j3, and close 
to the critical point one expects finite size corrections as 
long as L or f3 are smaller than the spatial correlation 
length or imaginary correlation time, respectively. We 
used (3 < 1000 and averaged over at least 512 samples. 

In Fig. [l] we show our results for d/z'(h) in the 
Griffiths-McCoy region. For ft->oowe expect d/z'(h) — > 
1, since this is the result for isolated spins in random fields 
with non- vanishing probability weight at hi = 0. The 
more interesting limit is /i — ^ /i c . The data are well QGp,- 
patible with limh—>h c z'(h) = oo, as in one-dini£iisionQ'll3. 
The average susceptibility [%] av diverges likqifl [x] av ~ 
T d/z'(h)-i for T _^ q Hencej if lim^^ z '(h) is univer- 
sal, i.e. does not depend on the details of the disorder, 
as is the case in 1-d, [x]av diverges quite generally in a 
range about the quantum critical point for systems with 
Ising symmetry. 

Next we describe our results for critical exponents us- 
ing the discrete imaginary time algorithm, for which we 
studied sizes up to L = 48 and L T = 2048. We found that 
no more than 100 sweeps were required for equilibration, 
even for the largest lattices. At least 1000 samples were 
averaged over. 

We locate the T = critical point hy a method al- 
ready used for the quantum spin glasstj. We compute 
the Binder ratio, Eq. (^) , which (assuming, for now, 
that z is finite) has the finite-size scaling form 

ffav =?(^L 1 /^ ! L T /L z ), (9) 

where S — T c i — T c c ; , with the value of T c i at criti- 
cality. For fixed L, g &v has a peak as a function of L T . 
At the critical point, Tf, the peak height is independent 
of L and the values of L T at the maximum, L™ ax , vary 
as L z . Furthermore, a plot of g av against L T /L" lax at 
the critical point, which has the advantage of not need- 
ing a value for z, should collapse the data. We see in 
Fig. H that this does not happen. Rather the curves 
clearly become broader for larger sizes. This is easy to 
understand since we know that for 1-d z is infinite and 
it is the log of the characteristic time which scales with 
a power of the length scale. This suggests that the scal- 
ing variable should be lnL r /lnL™ ax with lni™ ax ~ L^, 
say. A corresponding scaling plot is shown in the inset 
to Fig. ||. The data collapse for sizes L > 12, quite good 



for ip = 0.42, and not quite so good with the 1-d value, 
-0 = 1/2, though we would not claim that ip — 1/2 is 
ruled out. We have also performed analogous calcula- 
tions for one-dimensionE3 with very similar results. The 
close similarity of the data for 1-d and 2-d, and the broad- 
ening of the data in the main part of Fig. g, suggests that 
z is infinite in 2-d, as well as in 1-d. 
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FIG. 2. Results using the discrete imaginary time algo- 
rithm at the quantum critical point, T c i = Tf = 2.45. In the 
main figure the horizontal axis is L T /L™ ax where L™ ax is the 
value of L T at the peak. Note that the curves do not scale 
but rather get broader for larger sizes, indicating activated 
scaling, z = oo, The data for L = 48 is slightly high which 
may indicate that the true value of Tf is a little higher. In 
the inset, the data for L > 12 is seen to scale quite well with 
the same form In Lr/L^ known to be exact in 1-d, but the 
value of ip = 0.42 (shown) works better than the 1-d value of 
V> = 1/2. 

Next we consider the equal time correlations at the 
criticaL point. Fig. || shows data for the average and 
typicaB correlations for spins separated by r — (L/2,0). 
For each value of L, we took L T such that g av is close 
to the peakEll shown in Fig. ||. According to finite size 
scaling, the dependence on L for a finite system should 
be the same as the dependence on r in a bulk system. 

The data in Fig. |^ shows that the average correla- 
tion function falls off with a power law, with fj about 
1.95, while the typical value falls off faster than a 
power law, (because of the downward curvature) consis- 
tent with a stretched exponential behavior of the form 
cxp(— const. L a ), with a ~ 1/3. The statistical er- 
rors (shown) are generally smaller than the size of the 
points, so the downward curvature is statistically sig- 
nificant. This behavior is of the form expected in one 
dimensioncl, except that there a = a result which is 
reproduced by our 1-d simulationsEJ. Moreover, in one 
dimension iK= o~, and this relation also holds in higher 
dimensionstJ provided that the fixed point is similar, i.e. 
has infinitely strong disorder. Our data is compatible 
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FIG. 3. The main figure shows the average and typicalEl 
correlations between spins L/2 apart at the critical point. The 
average falls off with a power law, and a fit gives a slope of 
— fj with fj = 1.95. However, the curvature of the data for 
the typical correlation function shows that this falls off faster 
than a power law. The inset shows that the data is consistent 
with a stretched exponential form, exp(— const. L"), but with 
cr = l/3 rather than the value of 1/2 foundu in 1-d. 



with this result though neither a nor ip are determicnd 
with precision. Additional results, including the whole 
distribution of correlation functions, will be presented in 
a separate publication. 

To conclude, we have found a strong similarity between 
the critical behavior of the random transverse field Ising 
model in one and two dimensions. In particular, both 
z and lim/ l _ > ^ e z'(K) arc oafinitc. Previous simulations 
on quantum spin glassesll^l (for which the mean of the 
distribution of the Jij is zero) in two dimensions, found 
these quantities to be apparently finite. However, it is 
plausible that the asymptotic result should be infinite 
also for quantum spin glasses, and that the finite result 
found is an artifact of the smaller sizes used there. (Those 
studies also used a non-random transverse field.) 

Two of us (NK and HR) have implemented numerically 
for d = 2 a natural generalization of the renormalization 
group procedure used in Ref. || for d = 1 . After this work 
was completed, we heard that S.-C. Mau, O. Motrunich 
and D. A. Huse (private communication) used a simi- 
lar approach and observed a flow to the infinite disorder 
critical fixed point, just as in d = 1. 
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